
/***
 * Program to subtract background intensity from a time-density curve.
 * Usage: $ bgsub p.txt p.bgr p.dsa
 * p.txt: Original TDC.
 * p.bgr: Background ROI as "topleftX topleftY width height"
 * p.dsa: Image sequence.
 * 
 ***/

#include <iostream>
#include <vector>
using namespace std;

#include "dsaio.h"
 
int main ( int argc, char **argv )
{
  // process command line args
  if ( argc != 4 )
  {
    cout << "Usage: " << argv[ 0 ] << " p.txt p.bgr p.dsa" << endl;
	 return -1;
  }
  char *txtfile = argv[ 1 ];
  char *bgrfile = argv[ 2 ];
  char *dsafile = argv[ 3 ];
  
  // read input TDC
  // cerr << "  read input - ";
  vector< double > time_vector;
  vector< double > density_vector;
  ifstream txtinput;
  txtinput.open ( txtfile );
  double time, density;
  while ( txtinput >> time && txtinput >> density )
  {
    time_vector.push_back ( time );
	 density_vector.push_back ( density );
  }
  txtinput.close();
  // cerr << time_vector.size() << " data read" << endl;
  
  // read roi
  // cerr << "  read roi (";
  ifstream bgrinput;
  bgrinput.open ( bgrfile );
  int min_x, min_y, width, height;
  bgrinput >> min_x;
  bgrinput >> min_y;
  bgrinput >> width;
  bgrinput >> height;
  // cerr << min_x << ", " << min_y << ", "
  //     << width << ", " << height << ")" << endl;
  
  // compute and subtract background values
  // cerr << "  iterate frames" << endl;
  typedef itk::ImageRegionConstIterator< DsaFrame > DsaConstIterator;
  vector< double > bgs_vector;
  for ( int fr = 0; fr < time_vector.size(); ++ fr )
  {
    DsaFrame::Pointer frame = dsa_get_frame ( dsafile, fr );
    DsaFrame::RegionType roi;
    DsaFrame::RegionType::IndexType roiStart;
    DsaFrame::RegionType::SizeType roiSize;
    roiStart[ 0 ] = min_x;
    roiStart[ 1 ] = min_y;
    roiSize[ 0 ] = width;
    roiSize[ 1 ] = height;
    roi.SetSize( roiSize );
    roi.SetIndex( roiStart );
    DsaConstIterator it ( frame, roi );
    double avg = 0.0;
    double denom = 0.0;
    for ( it.GoToBegin(); ! it.IsAtEnd(); ++ it ) // iterate pixels
    {
      avg += (double)( it.Get() );
      denom += 1.0;
    }
    avg /= denom;
    bgs_vector.push_back ( density_vector[ fr ] - avg );
  }
  
  // cerr << time_vector.size() << " " << density_vector.size() << endl;
  
  // write output
  // put zero level back where it was in the original TDC
  // cerr << "  write output" << endl;
  double zero_level = density_vector[ 0 ] - bgs_vector[ 0 ];
  ofstream output;
  output.open ( txtfile );
  for ( int i = 0; i < time_vector.size(); ++ i )
  {
    output << time_vector[ i ] << " " << bgs_vector[ i ] + zero_level << endl;
  }
  output.close();
  
  return 0;
}
